package mirabilis.org.util;

import random.CommonState;
import mirabilis.org.IEvaluator;
import mirabilis.org.config.Configuration;
import mirabilis.org.simplices.Coordinates;
import mirabilis.org.simplices.Simplex;
import mirabilis.org.simplices.Vertex;

public class Initialization {
	
	public static Vertex[] init_simplex(IEvaluator evaluator){
		Vertex[] simplex = new Vertex[Configuration.N+1];
		
		double[][] csimplex = new double[Configuration.N+1][Configuration.N];
		
		//centroid
		for(int i=0;i<Configuration.N;i++){
			double lbound = Configuration.bounds[i*2];
			double ubound = Configuration.bounds[(i*2)+1];
			double center =  (lbound + ubound) / 2.0;
			
			//center = CommonState.r.nextDouble()*3.0;
			csimplex[0][i] = center;
		}
		
		//Random in innerbox
		
		
		
		simplex[0] = new Vertex(new Coordinates(csimplex[0]));
		simplex[0].setFitness(evaluator.evaluate(csimplex[0]));
		
		
		
		//Make all points initially equal to the first one
		for(int i=1;i<Configuration.N+1;i++){
			System.arraycopy(csimplex[0], 0, csimplex[i], 0, Configuration.N);
		}
		
		//Move in every orthogonal direction (except for the initial coordinate)
		// The movement is 1/4 of the bounded dimension either to left or right
		// we choose the best out of them.
		// Here we may be able to save one evaluation (we assume unimodality):
		// Consider the orthogonal line in which O is the center
		// L is the Left movement and R the Right one
		//     -----L-----O-----R-----
		// if f(R) is better than f(O), it is obvious that f(L) will be worse than f(R)
		// when that is the case we save the evaluation of f(L)
		for(int i=1;i<Configuration.N+1;i++){
			int dim = i-1;
			double lbound = Configuration.bounds[dim*2];
			double ubound = Configuration.bounds[(dim*2)+1];
			
			double move = (ubound-lbound)/4.0;
			csimplex[i][dim] = csimplex[0][dim] + move;
			Vertex movingRightVertex = new Vertex(new Coordinates(csimplex[i]));
			movingRightVertex.setFitness(evaluator.evaluate(csimplex[i]));

			if(!movingRightVertex.better(simplex[0])){
				csimplex[i][dim] = csimplex[0][dim] - move;
				Vertex movingLeftVertex = new Vertex(new Coordinates(csimplex[i]));
				movingLeftVertex.setFitness(evaluator.evaluate(csimplex[i]));
				if(movingLeftVertex.better(movingRightVertex)){
					simplex[i] = movingLeftVertex;
				}else{
					simplex[i] = movingRightVertex;
				}
			}else{
				simplex[i] = movingRightVertex;
			}		
		}
		
		return simplex;
		
	}
	
	/*
	  Purpose:

	    SIMPLEX_COORDINATES1 computes the Cartesian coordinates of simplex vertices.

	  Discussion:

	    The simplex will have its centroid at 0;

	    The sum of the vertices will be zero.

	    The distance of each vertex from the origin will be 1.

	    The length of each edge will be constant.

	    The dot product of the vectors defining any two vertices will be - 1 / N.
	    This also means the angle subtended by the vectors from the origin
	    to any two distinct vertices will be arccos ( - 1 / N ).

	  Licensing:

	    This code is distributed under the GNU LGPL license.

	  Modified:

	    20 September 2010

	  Author:

	    John Burkardt

	  Parameters:

	    Input, int N, the spatial dimension.

	    Output, double SIMPLEX_COORDINATES1[N*(N+1)], the coordinates of the vertices
	    of a simplex in N dimensions.  
	*/
	
	public static double [] simplex_coordinates ( int n )
	{
	  int i;
	  int ii;
	  int j;
	  double s;
	  double [] x;

	  x = r8mat_zero_new ( n, n + 1 );

	  for ( i = 0; i < n; i++ )
	  {
	/*
	  Set X(I,I) so that sum ( X(1:I,I)**2 ) = 1.
	*/
	    s = 0.0;
	    for ( ii = 0; ii < i; ii++ )
	    {
	      s = s + x[ii+i*n] * x[ii+i*n];
	    }
	    x[i+i*n] = Math.sqrt ( 1.0 - s );
	/*
	  Set X(I,J) for J = I+1 to N+1 by using the fact that XI dot XJ = - 1 / N 
	*/
	    for ( j = i + 1; j < n + 1; j++ )
	    {
	      s = 0.0;
	      for ( ii = 0; ii < i; ii++ )
	      {
	        s = s + x[ii+i*n] * x[ii+j*n];
	      }
	      x[i+j*n] = ( - 1.0 / ( double ) ( n ) - s ) / x[i+i*n];
	    }
	  }

	  return x;
	}
	
	
	public static double simplex_volume ( int n, double x[] )

	/******************************************************************************/
	/*
	  Purpose:

	    SIMPLEX_VOLUME computes the volume of a simplex.

	  Licensing:

	    This code is distributed under the GNU LGPL license.

	  Modified:

	    19 September 2010

	  Author:

	    John Burkardt

	  Parameters:

	    Input, int N, the spatial dimension.

	    Input, double X[N*(N+1)], the coordinates of the vertices
	    of a simplex in N dimensions.  

	    Output, double SIMPLEX_VOLUME, the volume of the simplex.
	*/
	{
	  double []a;
	  double det;
	  int i;
	  int j;
	  double volume;

	  a = new double[n*n];
	  
	  for ( j = 0; j < n; j++ )
	  {
	    for ( i = 0; i < n; i++ )
	    {
	      a[i+j*n] = x[i+j*n];
	    }
	  }
	  for ( j = 0; j < n; j++ )
	  {
	    for ( i = 0; i < n; i++ )
	    {
	      a[i+j*n] = a[i+j*n] - x[i+n*n];
	    }
	  }

	  det = r8mat_det ( n, a );

	  volume = r8_abs ( det );
	  for ( i = 1; i <= n; i++ )
	  {
	    volume = volume / ( double ) ( i );
	  }


	  return volume;
	}
	
	
	private static double [] r8mat_zero_new ( int m, int n )

	{
	  double [] a =new double[m*n];

	  for ( int j = 0; j < n; j++ )
	  {
	    for ( int i = 0; i < m; i++ )
	    {
	      a[i+j*m] = 0.0;
	    }
	  }
	  return a;
	}
	
	private static double r8mat_det ( int n, double a[] )

	/******************************************************************************/
	/*
	  Purpose:

	    R8MAT_DET computes the determinant of an R8MAT.

	  Discussion:

	    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector 
	    in column-major order.

	  Licensing:

	    This code is distributed under the GNU LGPL license. 

	  Modified:

	    14 May 2010

	  Author:

	    Original FORTRAN77 version by Helmut Spaeth.
	    C version by John Burkardt.
	    Java version by Juanlu Jimenez.

	  Reference:

	    Helmut Spaeth,
	    Cluster Analysis Algorithms
	    for Data Reduction and Classification of Objects,
	    Ellis Horwood, 1980, page 125-127.

	  Parameters:

	    Input, int N, the order of the matrix.

	    Input, double A[N*N], the matrix whose determinant is desired.

	    Output, double R8MAT_DET, the determinant of the matrix.
	*/
	{
	  double []b;
	  double det;
	  int i;
	  int j;
	  int k;
	  int kk;
	  int m;
	  double temp;

	  b = new double[n*n];

	  for ( j = 0; j < n; j++ )
	  {
	    for ( i = 0; i < n; i++ )
	    {
	      b[i+j*n] = a[i+j*n];
	    }
	  }

	  det = 1.0;

	  for ( k = 1; k <= n; k++ )
	  {
	    m = k;
	    for ( kk = k+1; kk <= n; kk++ )
	    {
	      if ( r8_abs ( b[m-1+(k-1)*n] ) < r8_abs ( b[kk-1+(k-1)*n] ) )
	      {
	        m = kk;
	      }
	    }

	    if ( m != k )
	    {
	      det = -det;

	      temp = b[m-1+(k-1)*n];
	      b[m-1+(k-1)*n] = b[k-1+(k-1)*n];
	      b[k-1+(k-1)*n] = temp;
	    }

	    det = det * b[k-1+(k-1)*n];

	    if ( b[k-1+(k-1)*n] != 0.0 )
	    {
	      for ( i = k+1; i <= n; i++ )
	      {
	        b[i-1+(k-1)*n] = -b[i-1+(k-1)*n] / b[k-1+(k-1)*n];
	      }

	      for ( j = k+1; j <= n; j++ )
	      {
	        if ( m != k )
	        {
	          temp = b[m-1+(j-1)*n];
	          b[m-1+(j-1)*n] = b[k-1+(j-1)*n];
	          b[k-1+(j-1)*n] = temp;
	        }
	        for ( i = k+1; i <= n; i++ )
	        {
	          b[i-1+(j-1)*n] = b[i-1+(j-1)*n] + b[i-1+(k-1)*n] * b[k-1+(j-1)*n];
	        }
	      }
	    }
	  }
	  return det;
	}
	
	private static double r8_abs ( double x )

	/******************************************************************************/
	/*
	  Purpose:

	    R8_ABS returns the absolute value of an R8.

	  Licensing:

	    This code is distributed under the GNU LGPL license. 

	  Modified:

	    07 May 2006

	  Author:

	    John Burkardt

	  Parameters:

	    Input, double X, the quantity whose absolute value is desired.

	    Output, double R8_ABS, the absolute value of X.
	*/
	{
	  double value;

	  if ( 0.0 <= x )
	  {
	    value = + x;
	  } 
	  else
	  {
	    value = - x;
	  }
	  return value;
	}
	
	
	public static double volume_regular(int n){
		
		double fact=1;
		for (int i = 2; i <= n; i++ )
			fact *= i;
		
		double side = Math.sqrt( (2.0*(n+1))/n );
		
		//return (Math.sqrt(n+1)/(fact*Math.sqrt(Math.pow(2.0, n))))*Math.pow(1.0/Math.sqrt(2), n+1);
		return (Math.sqrt(n+1)/(fact*Math.sqrt(Math.pow(2.0, n)))) * Math.pow(side, n);
		
	}
	
	public static double magnitude_of_scale(int n, double targeted_volume){
		
		double fact=1.0;
		for (int i = 2; i <= n; i++ )
			fact *= i;
		
		
		
		double o = Math.pow((targeted_volume * fact) / ( Math.sqrt((n+1.0)/Math.pow(2.0, n)) ),1.0/n) / Math.sqrt((2.0*(n+1))/ (n)); 
		
		return o;
	}
	
	public static double[] scale_regular_simplex(double o, double []x){
		double [] scaled = new double[x.length];
		
		for (int i=0;i<x.length;i++){
			scaled[i] = x[i]*o;
		}
		return scaled;
	}
	
	public static double volume_constrained_hypercube(int n,double[] lu_const)
	{
		double v = 1;
		//Every side of the hypercube is given as a pair of lower and upper constraints
		for(int i=0;i<2*n;i+=2){
			v *= (lu_const[i+1]-lu_const[i]);
		}
		return v;
	}
	
	//lu stands for lower and upper bounds
	public static double[] toroidal_coordinates(double[] coor, double[] lu_const){
		
		double[] tor_coor = new double[coor.length];
		
		// Uncomment for a real torus
		
//		for(int i=0;i<coor.length;i++){
//			double aux = ((coor[i]-lu_const[(i*2)])%(lu_const[(i*2)+1]-lu_const[(i*2)]));
//			if (aux > 0)
//				tor_coor[i] = aux+lu_const[(i*2)];
//			else
//				tor_coor[i] = lu_const[(i*2)+1] + aux;			
//		}

		// Stay in the bounds
		for(int i=0;i<coor.length;i++){

			if (coor[i] < lu_const[(i*2)])// smaller than lower bound
				tor_coor[i] = lu_const[i*2];
			else if(coor[i] > lu_const[(i*2)+1])// larger than upper bound
				tor_coor[i] = lu_const[(i*2)+1];
			else tor_coor[i] = coor[i];
		}

		
		return tor_coor;
	}
	
	
	
	
	
	public static void main(String[] args) {
		
		int n = 3;
		
		double[] x = Initialization.simplex_coordinates(n);
		 
		
		
		
		for(int i=0; i<x.length;i++){
			if(i%n==0){
				System.out.println();
			
			}
			System.out.print(x[i]+" ");
			
		}
		System.out.println();
		
		System.out.println("-----");
		System.out.println("Volume: "+simplex_volume(n, x));
		System.out.println("Volume regular: "+volume_regular(n));
		System.out.println("-----");
		
		
		double o = Initialization.magnitude_of_scale(n, 4000000.0);
		
		double[] sc = Initialization.scale_regular_simplex(o, x);
		
		double[][] nsimplex = new double[n+1][n];
		int point=-1;
		for(int i=0; i<sc.length;i++){
			if(i%n==0){
				System.out.println();
				point++;
			}
			System.out.print(sc[i]+" ");
			nsimplex[point][i%n] = sc[i];
		}
		point++;
		System.out.println();
		
		System.out.println("-----");
		System.out.println("Volume: "+simplex_volume(n, sc));
		System.out.println("Volume regular: "+volume_regular(n));
		System.out.println("-----");
		
		//Hypercube

		int nn=3;
		double []H = {-100,100,-100,100,-100,100};
		
		System.out.println("Volume HC: "+volume_constrained_hypercube(nn, H));
		
		for(int i=0; i<point;i++){
			double[] coord_torus = toroidal_coordinates(nsimplex[i], H);
			System.out.println();
			for(int j=0;j<coord_torus.length;j++)
				System.out.print(coord_torus[j]+" ");
		}
		
		
	}

}
